*!version 6.0  14Oct2014 

capture mata mata drop rdvce()
mata
real matrix rdvce(real matrix X, real matrix y, real matrix z, real scalar p, real scalar h, real scalar matches, string vce, string kernel)
{
n = length(X)
p1 = p+1
sigma = J(n,1, .)
if (vce=="resid") {
	for (k=1; k<=n; k++) {
		cutoff = X[k]
		W = kweight(X,cutoff,h,kernel)
		W_p = select(W, W:> 0)
		X_p = select(X, W:> 0)
		y_p = select(y, W:> 0)
		z_p = select(z, W:> 0)
		XX  = J(length(W_p),p1,.)
		for (j=1; j<=p1; j++) {
			XX[.,j] = (X_p:-cutoff):^(j-1)
		}
	m_p_y = invsym(cross(XX,W_p,XX))*cross(XX,W_p,y_p)
	m_p_z = invsym(cross(XX,W_p,XX))*cross(XX,W_p,z_p)
	sigma[k] = (y[k] - m_p_y[1])*(z[k] - m_p_z[1])
	}
}
else  {
v = w = 0
for (k=1; k<=n; k++) {
	x_abs = abs(X :- X[k,1])
	diffx = select(x_abs, x_abs:> 0)
	diffy = select(y, x_abs:> 0)
	diffz = select(z, x_abs:> 0)
	minindex(diffx, matches, v, w)
	Ji=length(v)
	y_match_avg = mean(diffy[v])
	z_match_avg = mean(diffz[v])
	sigma[k] = (Ji/(Ji+1))*(y[k] :- y_match_avg):*(z[k] :- z_match_avg)
	}
}
return(sigma)
}
mata mosave rdvce(), replace
end

capture mata mata drop kweight()
mata
real matrix kweight(real matrix X, real scalar c, real scalar h, string kernel)
{
u = (X:-c)/h
	if (kernel=="epanechnikov" | kernel=="epa") {
	w = (0.75:*(1:-u:^2):*(abs(u):<=1))/h
	}
	else if (kernel=="uniform" | kernel=="uni") {
	w = (0.5:*(abs(u):<=1))/h
	}
	else {
	w = ((1:-abs(u)):*(abs(u):<=1))/h
	}	
return(w)	
}
mata mosave kweight(), replace
end

capture mata mata drop regconst()
mata
real matrix regconst(real scalar d, real scalar h)
{
d2 = 2*d+1
d1 = d+1
mu = J(d2, 1, 0)
mu[1] = 1
XX = J(d1,d1,0)
for (j=2; j<=d2; j++) {
i = j-1
	if (mod(j,2)==1) {
		mu[j] = (1/(i+1))*(h/2)^i
	}
}
for (j=1; j<=d1; j++) {
	XX[j,.] = mu[j::j+d]'
}
invXX =invsym(XX)
return(invXX)
}
mata mosave regconst(), replace
end

capture mata mata drop cvplot()
mata
void cvplot(
 real colvector y,
 real colvector x,
 | string scalar opts)
{
	real scalar n, N, Y, X

	n = rows(y)
	if (rows(x)!=n) _error(3200)
	N = st_nobs()
	if (N<n) st_addobs(n-N)
	st_store((1,n), Y=st_addvar("double", st_tempname()), y)
	st_store((1,n), X=st_addvar("double", st_tempname()), x)
	stata("twoway scatter " + st_varname(Y) + " " +
	 st_varname(X) + ", " + opts)
	if (N<n) st_dropobsin((N+1,n))
	st_dropvar((Y,X))
}
mata mosave cvplot(), replace
end






